#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <errno.h>
#include <cinttypes>
#include <algorithm>
#include <tuple>
#include <omp.h>
#include <iostream>
#include "kthread.h"
#include "kvec.h"
#include "kalloc.h"
#include "sdust.h"
#include "mmpriv.h"
#include "bseq.h"
#include "khash.h"

struct mm_tbuf_s {
	void *km;
	int rep_len, frag_gap;
};

mm_tbuf_t *mm_tbuf_init(void)
{
	mm_tbuf_t *b;
	b = (mm_tbuf_t*)calloc(1, sizeof(mm_tbuf_t));
	if (!(mm_dbg_flag & 1)) b->km = km_init();
	return b;
}

void mm_tbuf_destroy(mm_tbuf_t *b)
{
	if (b == 0) return;
	km_destroy(b->km);
	free(b);
}

void *mm_tbuf_get_km(mm_tbuf_t *b)
{
	return b->km;
}

static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
{
	int n_dreg, j, k, u = 0;
	const uint64_t *dreg;
	sdust_buf_t *sdb;
	if (sdust_thres <= 0) return n;
	sdb = sdust_buf_init(km);
	dreg = sdust_core((const uint8_t*)seq, l_seq, sdust_thres, 64, &n_dreg, sdb);
	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
		int32_t s = qpos - (span - 1), e = s + span;
		while (u < n_dreg && (int32_t)dreg[u] <= s) ++u;
		if (u < n_dreg && (int32_t)(dreg[u]>>32) < e) {
			int v, l = 0;
			for (v = u; v < n_dreg && (int32_t)(dreg[v]>>32) < e; ++v) { // iterate over LCRs overlapping this minimizer
				int ss = s > (int32_t)(dreg[v]>>32)? s : dreg[v]>>32;
				int ee = e < (int32_t)dreg[v]? e : (uint32_t)dreg[v];
				l += ee - ss;
			}
			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
		} else a[k++] = a[j];
	}
	sdust_buf_destroy(sdb);
	return k; // the new size
}

static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char* const* seqs, mm128_v *mv)
{
	int i, n, sum = 0;
	mv->n = 0;
	for (i = n = 0; i < n_segs; ++i) {
		size_t j;

		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv, mi);

		for (j = n; j < mv->n; ++j)
			mv->a[j].y += sum << 1;
		if (opt->sdust_thres > 0) // mask low-complexity minimizers
			mv->n = n + mm_dust_minier(km, mv->n - n, mv->a + n, qlens[i], seqs[i], opt->sdust_thres);
		sum += qlens[i], n = mv->n;
	}
}

#include "ksort.h"
#define heap_lt(a, b) ((a).x > (b).x)
KSORT_INIT(heap, mm128_t, heap_lt)

	typedef struct {
		uint32_t n;
		uint32_t q_pos, q_span;
		uint32_t seg_id:31, is_tandem:1;
		const uint64_t *cr;
	} mm_match_t;

static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
{
	int rep_st = 0, rep_en = 0, n_m;
	size_t i;
	mm_match_t *m;
	*n_mini_pos = 0;
	*mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
	m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t));
	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
		const uint64_t *cr;
		mm128_t *p = &mv->a[i];
		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
		int t;
		cr = mm_idx_get(mi, p->x>>8, &t);
		if (t >= max_occ) {
			int en = (q_pos >> 1) + 1, st = en - q_span;
			if (st > rep_en) {
				*rep_len += rep_en - rep_st;
				rep_st = st, rep_en = en;
			} else rep_en = en;
		} else {
			mm_match_t *q = &m[n_m++];
			q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32;
			q->is_tandem = 0;
			if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1;
			if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1;
			*n_a += q->n;
			(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q_pos>>1;
		}
	}
	*rep_len += rep_en - rep_st;
	*_n_m = n_m;
	return m;
}

static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const char *qname, int qlen, const mm_idx_t *mi, int *is_self)
{
	*is_self = 0;
	if (qname && (flag & (MM_F_NO_DIAG|MM_F_NO_DUAL))) {
		const mm_idx_seq_t *s = &mi->seq[r>>32];
		int cmp;
		cmp = strcmp(qname, s->name);
		if ((flag&MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen) {
			if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagnonal anchors
			if ((r&1) == (q->q_pos&1)) *is_self = 1; // this flag is used to avoid spurious extension on self chain
		}
		if ((flag&MM_F_NO_DUAL) && cmp > 0) // all-vs-all mode: map once
			return 1;
	}
	if (flag & (MM_F_FOR_ONLY|MM_F_REV_ONLY)) {
		if ((r&1) == (q->q_pos&1)) { // forward strand
			if (flag & MM_F_REV_ONLY) return 1;
		} else {
			if (flag & MM_F_FOR_ONLY) return 1;
		}
	}
	return 0;
}

static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
		int *n_mini_pos, uint64_t **mini_pos)
{
	int i, n_m, heap_size = 0;
	int64_t j, n_for = 0, n_rev = 0;
	mm_match_t *m;
	mm128_t *a, *heap;

	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);

	heap = (mm128_t*)kmalloc(km, n_m * sizeof(mm128_t));
	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));

	for (i = 0, heap_size = 0; i < n_m; ++i) {
		if (m[i].n > 0) {
			heap[heap_size].x = m[i].cr[0];
			heap[heap_size].y = (uint64_t)i<<32;
			++heap_size;
		}
	}
	ks_heapmake_heap(heap_size, heap);
	while (heap_size > 0) {
		mm_match_t *q = &m[heap->y>>32];
		mm128_t *p;
		uint64_t r = heap->x;
		int32_t is_self, rpos = (uint32_t)r >> 1;
		if (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) {
			if ((r&1) == (q->q_pos&1)) { // forward strand
				p = &a[n_for++];
				p->x = (r&0xffffffff00000000ULL) | rpos;
				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
			} else { // reverse strand
				p = &a[(*n_a) - (++n_rev)];
				p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
			}
			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
			if (is_self) p->y |= MM_SEED_SELF;
		}
		// update the heap
		if ((uint32_t)heap->y < q->n - 1) {
			++heap[0].y;
			heap[0].x = m[heap[0].y>>32].cr[(uint32_t)heap[0].y];
		} else {
			heap[0] = heap[heap_size - 1];
			--heap_size;
		}
		ks_heapdown_heap(0, heap_size, heap);
	}
	kfree(km, m);
	kfree(km, heap);

	// reverse anchors on the reverse strand, as they are in the descending order
	for (j = 0; j < n_rev>>1; ++j) {
		mm128_t t = a[(*n_a) - 1 - j];
		a[(*n_a) - 1 - j] = a[(*n_a) - (n_rev - j)];
		a[(*n_a) - (n_rev - j)] = t;
	}
	if (*n_a > n_for + n_rev) {
		memmove(a + n_for, a + (*n_a) - n_rev, n_rev * sizeof(mm128_t));
		*n_a = n_for + n_rev;
	}
	return a;
}

static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
		int *n_mini_pos, uint64_t **mini_pos)
{
	int i, n_m;
	mm_match_t *m;
	mm128_t *a;
	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
	for (i = 0, *n_a = 0; i < n_m; ++i) {
		mm_match_t *q = &m[i];
		const uint64_t *r = q->cr;
		uint32_t k;
		for (k = 0; k < q->n; ++k) {
			int32_t is_self, rpos = (uint32_t)r[k] >> 1;
			mm128_t *p;
			if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue;
			p = &a[(*n_a)++];
			if ((r[k]&1) == (q->q_pos&1)) { // forward strand
				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
			} else { // reverse strand
				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
			}
			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
			if (is_self) p->y |= MM_SEED_SELF;
		}
	}
	kfree(km, m);
	radix_sort_128x(a, a + (*n_a));
	return a;
}

static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
{
	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
		else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
		if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains
			mm_join_long(km, opt, qlen, n_regs, regs, a);
	}
}

static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, int *n_regs, mm_reg1_t *regs, mm128_t *a)
{
	if (!(opt->flag & MM_F_CIGAR)) return regs;
	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
		mm_set_sam_pri(*n_regs, regs);
	}
	return regs;
}

void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
{
	int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
	int max_chain_gap_qry, max_chain_gap_ref, min_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR);
	uint32_t hash;
	int64_t n_a;
	uint64_t *u, *mini_pos;
	mm128_t *a;
	mm128_v mv = {0,0,0};
	mm_reg1_t *regs0;
	km_stat_t kmst;

	//TODO: generalize this to n_segs > 1
	assert (n_segs == 1);		//deal with long reads (or asm contigs) only
	mm128_t **collect_a;
	int64_t *collect_n_a;
	mm_tbuf_t *b_master = b; //buffer for main thread entering this function

	//stage1: Pre-compute confident read alignments of substrings of input read
	//define new set of options for first stage
	//generate many candidate alignments to improve mapq estimation
	mm_mapopt_t opt2 = *opt;
	mm_mapopt_t *opt_2 = &opt2;
	opt_2->best_n = std::max(5, opt_2->best_n); //set minimum

	int countStartingPositions = 1 + std::ceil(qlens[0] * 1.0 / opt_2->suffixSampleOffset);
	collect_a = (mm128_t**)kmalloc(b->km, countStartingPositions * sizeof(mm128_t*)); 
	collect_n_a = (int64_t *)kmalloc(b->km, countStartingPositions * sizeof(int64_t));
	memset(collect_n_a, 0, countStartingPositions * sizeof(int64_t));

	//create a boolean vector to indicate what portion of read were mapped using MCASs
	int8_t* seqMapped = (int8_t *)kmalloc(b->km, qlens[0] * sizeof(int8_t));
	memset(seqMapped, 0, qlens[0] * sizeof(int8_t));

	//check if SVaware mode enabled and query length is sufficient
	if (opt_2->SVaware && qlens[0] >= opt_2->SVawareMinReadLength)
	{
		//parallelize single read alignment further for better load balance
#pragma omp parallel num_threads(OMP_PER_READ_THREADS)
		{
			//make all these variables private to openmp thread by redefining them
			int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
			uint32_t hash;
			int64_t n_a;
			uint64_t *u, *mini_pos;
			mm128_t *a;
			mm128_v mv = {0,0,0};
			mm_reg1_t *regs0;
			mm_tbuf_t *b = (mm_tbuf_t*)calloc(1, sizeof(mm_tbuf_t)); //omp thread local buffer
			b->km = km_init();
			int* sub_qlens = (int *)kmalloc(b->km, 1 * sizeof(int));
			char **sub_seqs = (char **) kmalloc(b->km, 1 * sizeof(char*));
			sub_seqs[0] = (char *)kmalloc(b->km, qlens[0] * sizeof(char));

#pragma omp for schedule(dynamic)
			for (int sub_begin = 0; sub_begin < qlens[0] + opt_2->suffixSampleOffset - 1; sub_begin += opt_2->suffixSampleOffset)
			{
				int suffix_id = sub_begin / opt_2->suffixSampleOffset;	//id for this string end-point

				bool mappingFound = false;
				int max_mapq_currentPos = 0;
				if (sub_begin >= qlens[0]) sub_begin = qlens[0]-1; //for last iter
				assert (sub_begin >= 0 && sub_begin < qlens[0]);

				for (int sub_len = opt_2->minPrefixLength; sub_len <= opt_2->maxPrefixLength; sub_len *= opt_2->prefixIncrementFactor)
				{
					//consider 'sub_len' bases to the right
					if (sub_begin + sub_len <= qlens[0])	//check substring end boundary limit
					{
						mv = {0,0,0};
						sub_qlens[0] = sub_len;
						memcpy (sub_seqs[0], &(seqs[0][sub_begin]), sub_len);

						for (i = 0, qlen_sum = 0; i < n_segs; ++i)
							qlen_sum += sub_qlens[i], n_regs[i] = 0, regs[i] = 0, n_regs0 = 0;

						if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) break;
						if (opt_2->max_qlen > 0 && qlen_sum > opt_2->max_qlen) break;

						hash  = qname? __ac_X31_hash_string(qname) : 0;
						hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt_2->seed);
						hash  = __ac_Wang_hash(hash);

						collect_minimizers(b->km, opt_2, mi, n_segs, sub_qlens, sub_seqs, &mv);
						if (opt_2->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_2, opt_2->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
						else a = collect_seed_hits(b->km, opt_2, opt_2->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

						if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
							fprintf(stderr, "RS\t%d\n", rep_len);
							for (i = 0; i < n_a; ++i)
								fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
										i == 0? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
						}

						// set max chaining gap on the query and the reference sequence
						if (is_sr)
							max_chain_gap_qry = qlen_sum > opt_2->max_gap? qlen_sum : opt_2->max_gap;
						else max_chain_gap_qry = opt_2->max_gap;

						if (opt_2->max_gap_ref > 0) {
							max_chain_gap_ref = opt_2->max_gap_ref; // always honor mm_mapopt_2_t::max_gap_ref if set
						} else if (opt_2->max_frag_len > 0) {
							max_chain_gap_ref = opt_2->max_frag_len - qlen_sum;
							if (max_chain_gap_ref < opt_2->max_gap) max_chain_gap_ref = opt_2->max_gap;
						} else max_chain_gap_ref = opt_2->max_gap;

						if (opt_2->min_gap_ref < max_chain_gap_ref)
							min_chain_gap_ref = opt_2->min_gap_ref;
						else min_chain_gap_ref = max_chain_gap_ref;

						a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_2->bw, opt_2->max_chain_skip, opt_2->max_chain_iter, opt_2->min_cnt, opt_2->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);

						if (opt_2->max_occ > opt_2->mid_occ && rep_len > 0) {
							int rechain = 0;
							if (n_regs0 > 0) { // test if the best chain has all the segments
								int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
								for (i = 0; i < n_regs0; ++i) { // find the best chain
									if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
									off += (uint32_t)u[i];
								}
								for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
									if ((a[max_off+i].y&MM_SEED_SEG_MASK) != (a[max_off+i-1].y&MM_SEED_SEG_MASK))
										++n_chained_segs;
								if (n_chained_segs < n_segs)
									rechain = 1;
							} else rechain = 1;
							if (rechain) { // redo chaining with a higher max_occ threshold
								kfree(b->km, a);
								kfree(b->km, u);
								kfree(b->km, mini_pos);
								if (opt_2->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_2, opt_2->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
								else a = collect_seed_hits(b->km, opt_2, opt_2->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
								a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_2->bw, opt_2->max_chain_skip, opt_2->max_chain_iter, opt_2->min_cnt, opt_2->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
							}
						}
						b->frag_gap = max_chain_gap_ref;
						b->rep_len = rep_len;

						regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);

						if (mm_dbg_flag & MM_DBG_PRINT_SEED)
							for (j = 0; j < n_regs0; ++j)
								for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i)
									fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
											i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));

						chain_post(opt_2, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
						if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);

						if (n_segs == 1) { // uni-segment
							regs0 = align_regs(opt_2, mi, b->km, sub_qlens[0], sub_seqs[0], &n_regs0, regs0, a);
							mm_set_mapq(b->km, n_regs0, regs0, opt_2->min_chain_score, opt_2->a, rep_len, is_sr);
							n_regs[0] = n_regs0, regs[0] = regs0;
						} else { // multi-segment
							mm_seg_t *seg;
							seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
							free(regs0);
							for (i = 0; i < n_segs; ++i) {
								mm_set_parent(b->km, opt_2->mask_level, opt_2->mask_len, n_regs[i], regs[i], opt_2->a * 2 + opt_2->b, opt_2->flag&MM_F_HARD_MLEVEL, opt_2->alt_drop); // update mm_reg1_t::parent
								regs[i] = align_regs(opt_2, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
								mm_set_mapq(b->km, n_regs[i], regs[i], opt_2->min_chain_score, opt_2->a, rep_len, is_sr);
							}
							mm_seg_free(b->km, n_segs, seg);
							if (n_segs == 2 && opt_2->pe_ori >= 0 && (opt_2->flag&MM_F_CIGAR))
								mm_pair(b->km, max_chain_gap_ref, opt_2->pe_bonus, opt_2->a * 2 + opt_2->b, opt_2->a, qlens, n_regs, regs); // pairing
						}

						int mostPromisingMapping = -1;
						int max_mapq_fragment = 0;

						//For valid mapping, save anchors 
						for (j = 0; j < n_regs0; ++j)
						{
							max_mapq_fragment = std::max ((int32_t)regs0[j].mapq, max_mapq_fragment);
							max_mapq_currentPos = std::max (max_mapq_fragment, max_mapq_currentPos);

							//Check for high confidence (mapq), length
							if (regs0[j].mapq >= opt_2->min_mapq && regs0[j].blen >= opt_2->min_qcov * sub_len && regs0[j].cnt > 0)
							{
								mappingFound = true;
								mostPromisingMapping = j;
								collect_n_a[suffix_id] = regs0[j].cnt;

								if (mm_dbg_flag & MM_DBG_POLISH)
								{
									//print MCAS information in paf-like  format, helpful for debugging & dot-plotting MCAS alignments
									fprintf(stderr, "PO\t%s %d %d %d %c %s %d %d %d %d %d %d %d [FOUND] \n", qname, qlens[0], sub_begin + regs0[j].qs, sub_begin + regs0[j].qe, "+-"[regs0[j].rev] , mi->seq[regs0[j].rid].name, mi->seq[regs0[j].rid].len, regs0[j].rs, regs0[j].re, regs0[j].mapq, suffix_id, sub_begin, sub_len);
								}

								break;		
							}
						}

						if ((mm_dbg_flag & MM_DBG_POLISH) && !mappingFound)
							fprintf(stderr, "PO\tqname:%s, suffid:%d, begin:%d, len:%d, max_mapq:%d, n_regs0:%d [NONE FOUND] \n", qname, suffix_id, sub_begin, sub_len, max_mapq_fragment, n_regs0);

						if (mappingFound)
						{
							assert (collect_n_a[suffix_id] > 0);
							assert (mostPromisingMapping >= 0);

#pragma omp critical
							{
								collect_a[suffix_id] = (mm128_t*)kmalloc(b_master->km, collect_n_a[suffix_id] * sizeof(mm128_t));
							}
							j = mostPromisingMapping;

							for (i = 0; i < regs0[j].cnt; ++i)
							{
								mm128_t _a_ = a[i + regs0[j].as];

								//correct coordinates of each anchor while storing
								if (_a_.x >> 63) //reverse strand 
									_a_.y += (qlens[0] - sub_begin - sub_len);
								else
									_a_.y += sub_begin;
								collect_a[suffix_id][i] = _a_;					
							}

							//mark mapped interval in boolean vector
#pragma omp critical
							{
								for(i = sub_begin; i < sub_begin + sub_len; i++)
									seqMapped[i] = 1;
							}
						}

						for (j = 0; j < n_regs0; ++j) {free (regs0[j].p);}
						free (regs0);
						kfree(b->km, mv.a);
						kfree(b->km, a);
						kfree(b->km, u);
						kfree(b->km, mini_pos);

						if (mappingFound || !n_regs0)
							break;		// mappingFound-> found shortest prefix; !n_regs0-> no candidate
					}

					//consider 'sub_len' bases to the left
					if (sub_begin - sub_len + 1 >= 0)			//check substring start boundary limit
					{
						mv = {0,0,0};
						sub_qlens[0] = sub_len;
						memcpy (sub_seqs[0], &(seqs[0][sub_begin - sub_len +1]), sub_len);

						for (i = 0, qlen_sum = 0; i < n_segs; ++i)
							qlen_sum += sub_qlens[i], n_regs[i] = 0, regs[i] = 0, n_regs0 = 0;

						if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) break;
						if (opt_2->max_qlen > 0 && qlen_sum > opt_2->max_qlen) break;

						hash  = qname? __ac_X31_hash_string(qname) : 0;
						hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt_2->seed);
						hash  = __ac_Wang_hash(hash);

						collect_minimizers(b->km, opt_2, mi, n_segs, sub_qlens, sub_seqs, &mv);
						if (opt_2->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_2, opt_2->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
						else a = collect_seed_hits(b->km, opt_2, opt_2->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

						if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
							fprintf(stderr, "RS\t%d\n", rep_len);
							for (i = 0; i < n_a; ++i)
								fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
										i == 0? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
						}

						// set max chaining gap on the query and the reference sequence
						if (is_sr)
							max_chain_gap_qry = qlen_sum > opt_2->max_gap? qlen_sum : opt_2->max_gap;
						else max_chain_gap_qry = opt_2->max_gap;

						if (opt_2->max_gap_ref > 0) {
							max_chain_gap_ref = opt_2->max_gap_ref; // always honor mm_mapopt_2_t::max_gap_ref if set
						} else if (opt_2->max_frag_len > 0) {
							max_chain_gap_ref = opt_2->max_frag_len - qlen_sum;
							if (max_chain_gap_ref < opt_2->max_gap) max_chain_gap_ref = opt_2->max_gap;
						} else max_chain_gap_ref = opt_2->max_gap;

						if (opt_2->min_gap_ref < max_chain_gap_ref)
							min_chain_gap_ref = opt_2->min_gap_ref;
						else min_chain_gap_ref = max_chain_gap_ref;

						a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_2->bw, opt_2->max_chain_skip, opt_2->max_chain_iter, opt_2->min_cnt, opt_2->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);

						if (opt_2->max_occ > opt_2->mid_occ && rep_len > 0) {
							int rechain = 0;
							if (n_regs0 > 0) { // test if the best chain has all the segments
								int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
								for (i = 0; i < n_regs0; ++i) { // find the best chain
									if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
									off += (uint32_t)u[i];
								}
								for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
									if ((a[max_off+i].y&MM_SEED_SEG_MASK) != (a[max_off+i-1].y&MM_SEED_SEG_MASK))
										++n_chained_segs;
								if (n_chained_segs < n_segs)
									rechain = 1;
							} else rechain = 1;
							if (rechain) { // redo chaining with a higher max_occ threshold
								kfree(b->km, a);
								kfree(b->km, u);
								kfree(b->km, mini_pos);
								if (opt_2->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_2, opt_2->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
								else a = collect_seed_hits(b->km, opt_2, opt_2->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
								a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_2->bw, opt_2->max_chain_skip, opt_2->max_chain_iter, opt_2->min_cnt, opt_2->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
							}
						}
						b->frag_gap = max_chain_gap_ref;
						b->rep_len = rep_len;

						regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);

						if (mm_dbg_flag & MM_DBG_PRINT_SEED)
							for (j = 0; j < n_regs0; ++j)
								for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i)
									fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
											i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));

						chain_post(opt_2, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
						if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);

						if (n_segs == 1) { // uni-segment
							regs0 = align_regs(opt_2, mi, b->km, sub_qlens[0], sub_seqs[0], &n_regs0, regs0, a);
							mm_set_mapq(b->km, n_regs0, regs0, opt_2->min_chain_score, opt_2->a, rep_len, is_sr);
							n_regs[0] = n_regs0, regs[0] = regs0;
						} else { // multi-segment
							mm_seg_t *seg;
							seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
							free(regs0);
							for (i = 0; i < n_segs; ++i) {
								mm_set_parent(b->km, opt->mask_level, opt->mask_len, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent
								regs[i] = align_regs(opt_2, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
								mm_set_mapq(b->km, n_regs[i], regs[i], opt_2->min_chain_score, opt_2->a, rep_len, is_sr);
							}
							mm_seg_free(b->km, n_segs, seg);
							if (n_segs == 2 && opt_2->pe_ori >= 0 && (opt_2->flag&MM_F_CIGAR))
								mm_pair(b->km, max_chain_gap_ref, opt_2->pe_bonus, opt_2->a * 2 + opt_2->b, opt_2->a, qlens, n_regs, regs); // pairing
						}

						int mostPromisingMapping = -1;
						int max_mapq_fragment = 0;

						//For valid mapping, save anchors 
						for (j = 0; j < n_regs0; ++j)
						{
							max_mapq_fragment = std::max ((int32_t)regs0[j].mapq, max_mapq_fragment);
							max_mapq_currentPos = std::max (max_mapq_fragment, max_mapq_currentPos);

							//Check for high confidence (mapq), length
							if (regs0[j].mapq >= opt_2->min_mapq && regs0[j].blen >= opt_2->min_qcov * sub_len && regs0[j].cnt > 0)
							{
								mappingFound = true;
								mostPromisingMapping = j;
								collect_n_a[suffix_id] = regs0[j].cnt;

								if (mm_dbg_flag & MM_DBG_POLISH)
								{
									//print MCAS information in paf-like  format, helpful for debugging & dot-plotting MCAS alignments
									fprintf(stderr, "PO\t%s %d %d %d %c %s %d %d %d %d %d %d %d [FOUND] \n", qname, qlens[0], sub_begin - sub_len + regs0[j].qs, sub_begin - sub_len + regs0[j].qe, "+-"[regs0[j].rev] , mi->seq[regs0[j].rid].name, mi->seq[regs0[j].rid].len, regs0[j].rs, regs0[j].re, regs0[j].mapq, suffix_id, sub_begin, -1 * sub_len);
								}

								break;	
							}
						}

						if ((mm_dbg_flag & MM_DBG_POLISH) && !mappingFound)
							fprintf(stderr, "PO\tqname:%s, suffid:%d, begin:%d, len:%d, max_mapq:%d, n_regs0:%d [NONE FOUND] \n", qname, suffix_id, sub_begin, -1 * sub_len, max_mapq_fragment, n_regs0);

						if (mappingFound)
						{
							assert (collect_n_a[suffix_id] > 0);
							assert (mostPromisingMapping >= 0);

#pragma omp critical
							{
								collect_a[suffix_id] = (mm128_t*)kmalloc(b_master->km, collect_n_a[suffix_id] * sizeof(mm128_t));
							}
							j = mostPromisingMapping;

							for (i = 0; i < regs0[j].cnt; ++i)
							{
								mm128_t _a_ = a[i + regs0[j].as];

								//correct coordinates of each anchor while storing
								if (_a_.x >> 63) //reverse strand 
									_a_.y += (qlens[0]-1) - sub_begin;  
								else
									_a_.y += sub_begin - sub_len + 1;		//offset of first base of substring
								collect_a[suffix_id][i] = _a_;					
							}

							//mark mapped interval in boolean vector
#pragma omp critical
							{
								for(i = sub_begin - sub_len +1; i <= sub_begin; i++)
									seqMapped[i] = 1;
							}
						}

						for (j = 0; j < n_regs0; ++j) {free (regs0[j].p);}
						free (regs0);
						kfree(b->km, mv.a);
						kfree(b->km, a);
						kfree(b->km, u);
						kfree(b->km, mini_pos);

						if (mappingFound || !n_regs0)
							break;		// mappingFound-> found shortest prefix; !n_regs0-> no candidate
					}
				}

				if ((mm_dbg_flag & MM_DBG_POLISH) && !mappingFound)
					fprintf(stderr, "PO\tqname:%s, begin:%d, max_mapq_currentPos:%d [NONE FOUND] \n", qname, sub_begin, max_mapq_currentPos);
			}

			//free openmp thread specific memory
			kfree(b->km, sub_qlens);
			kfree(b->km, sub_seqs[0]);
			kfree(b->km, sub_seqs);
			mm_tbuf_destroy(b);
		}
	}

	if (mm_dbg_flag & MM_DBG_POLISH)
	{
		int mappedcnt = 0;
		for (i = 0; i < qlens[0]; i++) if (seqMapped[i]) mappedcnt++;
		fprintf(stderr, "PO\tqname:%s, count of mapped query bases = %d among %d\n", qname, mappedcnt, qlens[0]);
	}

	//define new set of options for next stage
	//we can make stage 2 as sensitive as possible with very few seeds remaining
	mm_mapopt_t opt3 = *opt;
	mm_mapopt_t *opt_3 = &opt3;
	opt_3->zdrop_inv = std::min (opt->zdrop_inv, opt->stage2_zdrop_inv);
	opt_3->bw= std::max(opt->bw, opt->stage2_bw);

	//increased gap helps compensate for sometimes missing seeds along correct alignments
	opt_3->max_gap = std::max(opt->max_gap, opt->stage2_max_gap);

	//Re-run mapping with the above selected anchors
	{
		for (i = 0, qlen_sum = 0; i < n_segs; ++i)
			qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0, n_regs0 = 0;

		if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
		if (opt_3->max_qlen > 0 && qlen_sum > opt_3->max_qlen) return;

		hash  = qname? __ac_X31_hash_string(qname) : 0;
		hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt_3->seed);
		hash  = __ac_Wang_hash(hash);

		//Use anchors from our own analysis
		n_a = 0;
		for (i = 0; i < countStartingPositions; i++)
			n_a += collect_n_a[i];

		if ((mm_dbg_flag & MM_DBG_POLISH) && opt->SVaware)
			fprintf(stderr, "PO\tqname:%s, n_a (before filtering and checking for duplicates) :%" PRId64 "\n", qname, n_a);

		if (n_a)
		{
			//allocate sufficient memory
			a = (mm128_t*)kmalloc(b->km, n_a * sizeof(mm128_t));

			//set values of anchors
			int64_t n_a_counter = 0;
			for (i = 0; i < countStartingPositions; i++)
				for (j=0; j<collect_n_a[i]; j++)
					a[n_a_counter++] = collect_a[i][j];		

			//discard duplicate entries
			int64_t n_a_unique = 0;
			std::sort(a, a + n_a, [](const mm128_t &a, const mm128_t &b){return std::tie(a.x, a.y) < std::tie(b.x, b.y);});

			//traverse through the array elements
			for (i = 0; i < n_a;)
			{
				j = i;

				while (j < n_a && std::tie(a[i].x, a[i].y) == std::tie(a[j].x, a[j].y))
					j++;

				//j will increment at least by one here

				a[n_a_unique++] = a[i];
				i = j;
			}

			n_a = n_a_unique;

			if (mm_dbg_flag & MM_DBG_POLISH)
				fprintf(stderr, "PO\tqname:%s, n_a (after filtering and checking for duplicates) :%" PRId64 ", min_cnt:%d\n", qname, n_a, opt_3->min_cnt);

			//sort anchors by reference position before moving on
			radix_sort_128x(a, a + n_a);

			if (n_a < opt_3->min_cnt)  //insufficient no. of seeds
			{
				n_a = 0;	//reset to 0
				kfree(b->km, a);
			}
		}

	}

	//collect additional anchors from unmapped intervals
	{
		//if we have found MCAS-based anchors, but with a few unmapped read intervals
		int unmappedcnt = 0;
		for (i = 0; i < qlens[0]; i++) if (seqMapped[i]==0) unmappedcnt++;
		if (n_a > 0 && unmappedcnt > 0)
		{
			char **unmapped_seqs = (char **) kmalloc(b->km, 1 * sizeof(char*));
			unmapped_seqs[0] = (char *)kmalloc(b->km, qlens[0] * sizeof(char));
			for (i = 0; i < qlens[0]; i++)
			{
				if (seqMapped[i] > 0)
					unmapped_seqs[0][i] = 'N';
				else
					unmapped_seqs[0][i] = seqs[0][i];
			}

			if (mm_dbg_flag & MM_DBG_POLISH)
				fprintf(stderr, "PO\tqname:%s, n_a (before mapping the unmapped read substrings) :%" PRId64 "\n", qname, n_a);

			mm128_t *a_remaining;
			int64_t n_a_remaining;
			mv = {0,0,0};
			collect_minimizers(b->km, opt_3, mi, n_segs, qlens, unmapped_seqs, &mv);

			if (opt_3->flag & MM_F_HEAP_SORT)
				a_remaining = collect_seed_hits_heap(b->km, opt_3, opt_3->mid_occ, mi, qname, &mv, qlen_sum, &n_a_remaining, &rep_len, &n_mini_pos, &mini_pos);
			else
				a_remaining = collect_seed_hits(b->km, opt_3, opt_3->mid_occ, mi, qname, &mv, qlen_sum, &n_a_remaining, &rep_len, &n_mini_pos, &mini_pos);

			kfree(b->km, mv.a);
			kfree(b->km, mini_pos);

			int64_t n_a_whole = n_a_remaining + n_a;
			mm128_t *a_whole = (mm128_t*)kmalloc(b->km, n_a_whole * sizeof(mm128_t));

			for (i=0; i<n_a; i++)
			{
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
				a_whole[i] = a[i];
			}

			for (i=n_a; i<n_a_whole; i++)
			{
				a_whole[i] = a_remaining[i-n_a];
			}

			//sort anchors by reference position before moving on
			radix_sort_128x(a_whole, a_whole + n_a_whole);

			kfree(b->km, a);
			kfree(b->km, a_remaining);
			a = a_whole;
			n_a = n_a_whole;

			kfree(b->km, unmapped_seqs[0]);
			kfree(b->km, unmapped_seqs);

			if (mm_dbg_flag & MM_DBG_POLISH)
				fprintf(stderr, "PO\tqname:%s, n_a (after mapping the unmapped read substrings) :%" PRId64 "\n", qname, n_a);
		}
	}

	{
		if (!n_a) //MCAS-method couldn't be used
		{
			//go with the default route
			if ((mm_dbg_flag & MM_DBG_POLISH) && opt->SVaware)
				fprintf(stderr, "PO\tfalling back to default mapping algorithm for read: %s\n", qname);

			//revert to original parameters
			*opt_3 = *opt;

			mv = {0,0,0};
			collect_minimizers(b->km, opt_3, mi, n_segs, qlens, seqs, &mv);
			if (opt_3->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_3, opt_3->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
			else a = collect_seed_hits(b->km, opt_3, opt_3->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

			kfree(b->km, mv.a);
			kfree(b->km, mini_pos);
		}

		if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
			fprintf(stderr, "RS\t%d\n", rep_len);
			for (i = 0; i < n_a; ++i)
				fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
						i == 0? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
		}

		// set max chaining gap on the query and the reference sequence
		if (is_sr)
			max_chain_gap_qry = qlen_sum > opt_3->max_gap? qlen_sum : opt_3->max_gap;
		else max_chain_gap_qry = opt_3->max_gap;

		if (opt_3->max_gap_ref > 0) {
			max_chain_gap_ref = opt_3->max_gap_ref; // always honor mm_mapopt_3_t::max_gap_ref if set
		} else if (opt_3->max_frag_len > 0) {
			max_chain_gap_ref = opt_3->max_frag_len - qlen_sum;
			if (max_chain_gap_ref < opt_3->max_gap) max_chain_gap_ref = opt_3->max_gap;
		} else max_chain_gap_ref = opt_3->max_gap;

		if (opt_3->min_gap_ref < max_chain_gap_ref)
			min_chain_gap_ref = opt_3->min_gap_ref;
		else min_chain_gap_ref = max_chain_gap_ref;

		a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_3->bw, opt_3->max_chain_skip, opt_3->max_chain_iter, opt_3->min_cnt, opt_3->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);

		if (opt_3->max_occ > opt_3->mid_occ && rep_len > 0) {
			int rechain = 0;
			if (n_regs0 > 0) { // test if the best chain has all the segments
				int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
				for (i = 0; i < n_regs0; ++i) { // find the best chain
					if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
					off += (uint32_t)u[i];
				}
				for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
					if ((a[max_off+i].y&MM_SEED_SEG_MASK) != (a[max_off+i-1].y&MM_SEED_SEG_MASK))
						++n_chained_segs;
				if (n_chained_segs < n_segs)
					rechain = 1;
			} 
			else rechain = 1;
			if (rechain) { // redo chaining with a higher max_occ threshold
				kfree(b->km, a);
				kfree(b->km, u);
				//kfree(b->km, mini_pos); //already freed above
				if (opt_3->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt_3, opt_3->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
				else a = collect_seed_hits(b->km, opt_3, opt_3->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
				a = mm_chain_dp(max_chain_gap_ref, min_chain_gap_ref, max_chain_gap_qry, opt_3->bw, opt_3->max_chain_skip, opt_3->max_chain_iter, opt_3->min_cnt, opt_3->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
			}
		}
		b->frag_gap = max_chain_gap_ref;
		b->rep_len = rep_len;

		regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);

		if (mm_dbg_flag & MM_DBG_PRINT_SEED)
			for (j = 0; j < n_regs0; ++j)
				for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i)
					fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
							i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));

		chain_post(opt_3, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
		//This function generates lot of warnings
		/*if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);*/

		if (n_segs == 1) { // uni-segment
			regs0 = align_regs(opt_3, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
			mm_set_mapq(b->km, n_regs0, regs0, opt_3->min_chain_score, opt_3->a, rep_len, is_sr);
			n_regs[0] = n_regs0, regs[0] = regs0;
			//TODO: find a better way to compute mapping quality
		} else { // multi-segment
			mm_seg_t *seg;
			seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
			free(regs0);
			for (i = 0; i < n_segs; ++i) {
				mm_set_parent(b->km, opt_3->mask_level, opt_3->mask_len, n_regs[i], regs[i], opt_3->a * 2 + opt_3->b, opt_3->flag&MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent
				regs[i] = align_regs(opt_3, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
				mm_set_mapq(b->km, n_regs[i], regs[i], opt_3->min_chain_score, opt_3->a, rep_len, is_sr);
			}
			mm_seg_free(b->km, n_segs, seg);
			if (n_segs == 2 && opt_3->pe_ori >= 0 && (opt_3->flag&MM_F_CIGAR))
				mm_pair(b->km, max_chain_gap_ref, opt_3->pe_bonus, opt_3->a * 2 + opt_3->b, opt_3->a, qlens, n_regs, regs); // pairing
		}

		kfree(b->km, a);
		kfree(b->km, u);
		/*kfree(b->km, mini_pos);*/
		/*kfree(b->km, mv.a);*/
	}

	for (i = 0; i < countStartingPositions; i++)
		if (collect_n_a[i] > 0)
			kfree(b->km, collect_a[i]);

	kfree(b->km, collect_a);
	kfree(b->km, collect_n_a);
	kfree(b->km, seqMapped);

	if (b->km) {
		km_stat(b->km, &kmst);
		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
		if (kmst.largest > 1U<<28) {
			km_destroy(b->km);
			b->km = km_init();
		}
	}
}

mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
{
	mm_reg1_t *regs;
	mm_map_frag(mi, 1, &qlen, &seq, n_regs, &regs, b, opt, qname);
	return regs;
}

/**************************
 * Multi-threaded mapping *
 **************************/

typedef struct {
	int mini_batch_size, n_processed, n_threads, n_fp;
	const mm_mapopt_t *opt;
	mm_bseq_file_t **fp;
	const mm_idx_t *mi;
	kstring_t str;

	int n_parts;
	uint32_t *rid_shift;
	FILE *fp_split, **fp_parts;
} pipeline_t;

typedef struct {
	const pipeline_t *p;
	int n_seq, n_frag;
	mm_bseq1_t *seq;
	int *n_reg, *seg_off, *n_seg, *rep_len, *frag_gap;
	mm_reg1_t **reg;
	mm_tbuf_t **buf;
} step_t;

static void worker_for(void *_data, long i, int tid) // kt_for() callback
{
	step_t *s = (step_t*)_data;
	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
	const char *qseqs[MM_MAX_SEG];
	mm_tbuf_t *b = s->buf[tid];
	assert(s->n_seg[i] <= MM_MAX_SEG);
	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
	for (j = 0; j < s->n_seg[i]; ++j) {
		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
			mm_revcomp_bseq(&s->seq[off + j]);
		qlens[j] = s->seq[off + j].l_seq;
		qseqs[j] = s->seq[off + j].seq;
	}
	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
		for (j = 0; j < s->n_seg[i]; ++j) {
			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
			s->rep_len[off + j] = b->rep_len;
			s->frag_gap[off + j] = b->frag_gap;
		}
	} else {
		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
		for (j = 0; j < s->n_seg[i]; ++j) {
			s->rep_len[off + j] = b->rep_len;
			s->frag_gap[off + j] = b->frag_gap;
		}
	}
	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
			int k, t;
			mm_revcomp_bseq(&s->seq[off + j]);
			for (k = 0; k < s->n_reg[off + j]; ++k) {
				mm_reg1_t *r = &s->reg[off + j][k];
				t = r->qs;
				r->qs = qlens[j] - r->qe;
				r->qe = qlens[j] - t;
				r->rev = !r->rev;
			}
		}
}

static void merge_hits(step_t *s)
{
	int f, i, k0, k, max_seg = 0, *n_reg_part, *rep_len_part, *frag_gap_part, *qlens;
	void *km;
	FILE **fp = s->p->fp_parts;
	const mm_mapopt_t *opt = s->p->opt;

	km = km_init();
	for (f = 0; f < s->n_frag; ++f)
		max_seg = max_seg > s->n_seg[f]? max_seg : s->n_seg[f];
	qlens = CALLOC(int, max_seg + s->p->n_parts * 3);
	n_reg_part = qlens + max_seg;
	rep_len_part = n_reg_part + s->p->n_parts;
	frag_gap_part = rep_len_part + s->p->n_parts;
	for (f = 0, k = k0 = 0; f < s->n_frag; ++f) {
		k0 = k;
		for (i = 0; i < s->n_seg[f]; ++i, ++k) {
			int j, l, t, rep_len = 0;
			qlens[i] = s->seq[k].l_seq;
			for (j = 0, s->n_reg[k] = 0; j < s->p->n_parts; ++j) {
				mm_err_fread(&n_reg_part[j],    sizeof(int), 1, fp[j]);
				mm_err_fread(&rep_len_part[j],  sizeof(int), 1, fp[j]);
				mm_err_fread(&frag_gap_part[j], sizeof(int), 1, fp[j]);
				s->n_reg[k] += n_reg_part[j];
				if (rep_len < rep_len_part[j])
					rep_len = rep_len_part[j];
			}
			s->reg[k] = CALLOC(mm_reg1_t, s->n_reg[k]);
			for (j = 0, l = 0; j < s->p->n_parts; ++j) {
				for (t = 0; t < n_reg_part[j]; ++t, ++l) {
					mm_reg1_t *r = &s->reg[k][l];
					uint32_t capacity;
					mm_err_fread(r, sizeof(mm_reg1_t), 1, fp[j]);
					r->rid += s->p->rid_shift[j];
					if (opt->flag & MM_F_CIGAR) {
						mm_err_fread(&capacity, 4, 1, fp[j]);
						r->p = (mm_extra_t*)calloc(capacity, 4);
						r->p->capacity = capacity;
						mm_err_fread(r->p, r->p->capacity, 4, fp[j]);
					}
				}
			}
			mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
			mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
			if (!(opt->flag & MM_F_ALL_CHAINS)) {
				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]);
				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
			}
			mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));
		}
		if (s->n_seg[f] == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
			mm_pair(km, frag_gap_part[0], opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, &s->n_reg[k0], &s->reg[k0]);
	}
	free(qlens);
	km_destroy(km);
}

static void *worker_pipeline(void *shared, int step, void *in)
{
	int i, j, k;
	pipeline_t *p = (pipeline_t*)shared;
	if (step == 0) { // step 0: read sequences
		int with_qual = (!!(p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_NO_QUAL));
		int with_comment = !!(p->opt->flag & MM_F_COPY_COMMENT);
		int frag_mode = (p->n_fp > 1 || !!(p->opt->flag & MM_F_FRAG_MODE));
		step_t *s;
		s = (step_t*)calloc(1, sizeof(step_t));
		if (p->n_fp > 1) s->seq = mm_bseq_read_frag2(p->n_fp, p->fp, p->mini_batch_size, with_qual, with_comment, &s->n_seq);
		else s->seq = mm_bseq_read3(p->fp[0], p->mini_batch_size, with_qual, with_comment, frag_mode, &s->n_seq);
		if (s->seq) {
			s->p = p;
			for (i = 0; i < s->n_seq; ++i)
				s->seq[i].rid = p->n_processed++;

			//reshuffle based on length here, longer read first
			//NOTE: this would affect the ordering of reads in output
			{
				mm_bseq1_t *seq_copy = (mm_bseq1_t*) kmalloc(0, sizeof(mm_bseq1_t) * s->n_seq);
				std::vector< std::pair<int, int> > lengths;
				for (i = 0; i < s->n_seq; ++i)
					lengths.emplace_back (s->seq[i].l_seq, i);
				std::sort (lengths.begin(), lengths.end(), std::greater<std::pair<int,int>>());
				for (i = 0; i < s->n_seq; ++i) {
					int prev_id = lengths[i].second; //copy all pointers
					seq_copy[i].l_seq = s->seq[prev_id].l_seq;
					seq_copy[i].rid = s->seq[prev_id].rid;
					seq_copy[i].name = s->seq[prev_id].name;
					seq_copy[i].seq = s->seq[prev_id].seq;
					seq_copy[i].qual = s->seq[prev_id].qual;
					seq_copy[i].comment = s->seq[prev_id].comment;
				}
				free(s->seq);
				s->seq = seq_copy;
			}

			s->buf = (mm_tbuf_t**)calloc(p->n_threads, sizeof(mm_tbuf_t*));
			for (i = 0; i < p->n_threads; ++i)
				s->buf[i] = mm_tbuf_init();
			s->n_reg = (int*)calloc(5 * s->n_seq, sizeof(int));
			s->seg_off = s->n_reg + s->n_seq; // seg_off, n_seg, rep_len and frag_gap are allocated together with n_reg
			s->n_seg = s->seg_off + s->n_seq;
			s->rep_len = s->n_seg + s->n_seq;
			s->frag_gap = s->rep_len + s->n_seq;
			s->reg = (mm_reg1_t**)calloc(s->n_seq, sizeof(mm_reg1_t*));
			for (i = 1, j = 0; i <= s->n_seq; ++i)
				if (i == s->n_seq || !frag_mode || !mm_qname_same(s->seq[i-1].name, s->seq[i].name)) {
					s->n_seg[s->n_frag] = i - j;
					s->seg_off[s->n_frag++] = j;
					j = i;
				}
			return s;
		} else free(s);
	} else if (step == 1) { // step 1: map
		if (p->n_parts > 0) merge_hits((step_t*)in);
		else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
		return in;
	} else if (step == 2) { // step 2: output
		void *km = 0;
		step_t *s = (step_t*)in;
		const mm_idx_t *mi = p->mi;
		for (i = 0; i < p->n_threads; ++i) mm_tbuf_destroy(s->buf[i]);
		free(s->buf);
		if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init();
		for (k = 0; k < s->n_frag; ++k) {
			int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
			for (i = seg_st; i < seg_en; ++i) {
				mm_bseq1_t *t = &s->seq[i];
				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
					mm_err_fwrite(&s->n_reg[i],    sizeof(int), 1, p->fp_split);
					mm_err_fwrite(&s->rep_len[i],  sizeof(int), 1, p->fp_split);
					mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
					for (j = 0; j < s->n_reg[i]; ++j) {
						mm_reg1_t *r = &s->reg[i][j];
						mm_err_fwrite(r, sizeof(mm_reg1_t), 1, p->fp_split);
						if (p->opt->flag & MM_F_CIGAR) {
							mm_err_fwrite(&r->p->capacity, 4, 1, p->fp_split);
							mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
						}
					}
				} else if (s->n_reg[i] > 0) { // the query has at least one hit
					for (j = 0; j < s->n_reg[i]; ++j) {
						mm_reg1_t *r = &s->reg[i][j];
						assert(!r->sam_pri || r->id == r->parent);
						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
							continue;
						if (p->opt->flag & MM_F_OUT_SAM)
							mm_write_sam3(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
						else
							mm_write_paf3(&p->str, mi, t, r, km, p->opt->flag, s->rep_len[i]);
						mm_err_puts(p->str.s);
					}
				} else if ((p->opt->flag & MM_F_PAF_NO_HIT) || ((p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_SAM_HIT_ONLY))) { // output an empty hit, if requested
					if (p->opt->flag & MM_F_OUT_SAM)
						mm_write_sam3(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
					else
						mm_write_paf3(&p->str, mi, t, 0, 0, p->opt->flag, s->rep_len[i]);
					mm_err_puts(p->str.s);
				}
			}
			for (i = seg_st; i < seg_en; ++i) {
				for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p);
				free(s->reg[i]);
				free(s->seq[i].seq); free(s->seq[i].name);
				if (s->seq[i].qual) free(s->seq[i].qual);
				if (s->seq[i].comment) free(s->seq[i].comment);
			}
		}
		free(s->reg); free(s->n_reg); free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
		km_destroy(km);
		if (mm_verbose >= 3)
			fprintf(stderr, "[M::%s::%.3f*%.2f] mapped %d sequences\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), s->n_seq);
		free(s);
	}
	return 0;
}

static mm_bseq_file_t **open_bseqs(int n, const char **fn)
{
	mm_bseq_file_t **fp;
	int i, j;
	fp = (mm_bseq_file_t**)calloc(n, sizeof(mm_bseq_file_t*));
	for (i = 0; i < n; ++i) {
		if ((fp[i] = mm_bseq_open(fn[i])) == 0) {
			if (mm_verbose >= 1)
				fprintf(stderr, "ERROR: failed to open file '%s': %s\n", fn[i], strerror(errno));
			for (j = 0; j < i; ++j)
				mm_bseq_close(fp[j]);
			free(fp);
			return 0;
		}
	}
	return fp;
}

int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads)
{
	int i, pl_threads;
	pipeline_t pl;
	if (n_segs < 1) return -1;
	memset(&pl, 0, sizeof(pipeline_t));
	pl.n_fp = n_segs;
	pl.fp = open_bseqs(pl.n_fp, fn);
	if (pl.fp == 0) return -1;
	pl.opt = opt, pl.mi = idx;
	pl.n_threads = n_threads > 1? n_threads : 1;
	pl.mini_batch_size = opt->mini_batch_size;
	if (opt->split_prefix)
		pl.fp_split = mm_split_init(opt->split_prefix, idx);
	pl_threads = n_threads == 1? 1 : (opt->flag&MM_F_2_IO_THREADS)? 3 : 2;
	pl_threads = 1; 
	//TODO: this change helped avoid seg-faults on Phoenix cluster (figure out why)
	//GDB was indicating seg-faults in bseq.c 

	kt_pipeline(pl_threads, worker_pipeline, &pl, 3);

	free(pl.str.s);
	if (pl.fp_split) fclose(pl.fp_split);
	for (i = 0; i < pl.n_fp; ++i)
		mm_bseq_close(pl.fp[i]);
	free(pl.fp);
	return 0;
}

int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads)
{
	return mm_map_file_frag(idx, 1, &fn, opt, n_threads);
}

int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx)
{
	int i;
	pipeline_t pl;
	mm_idx_t *mi;
	if (n_segs < 1 || n_split_idx < 1) return -1;
	memset(&pl, 0, sizeof(pipeline_t));
	pl.n_fp = n_segs;
	pl.fp = open_bseqs(pl.n_fp, fn);
	if (pl.fp == 0) return -1;
	pl.opt = opt;
	pl.mini_batch_size = opt->mini_batch_size;

	pl.n_parts = n_split_idx;
	pl.fp_parts  = CALLOC(FILE*, pl.n_parts);
	pl.rid_shift = CALLOC(uint32_t, pl.n_parts);
	pl.mi = mi = mm_split_merge_prep(opt->split_prefix, n_split_idx, pl.fp_parts, pl.rid_shift);
	if (pl.mi == 0) {
		free(pl.fp_parts);
		free(pl.rid_shift);
		return -1;
	}
	for (i = n_split_idx - 1; i > 0; --i)
		pl.rid_shift[i] = pl.rid_shift[i - 1];
	for (pl.rid_shift[0] = 0, i = 1; i < n_split_idx; ++i)
		pl.rid_shift[i] += pl.rid_shift[i - 1];
	if (opt->flag & MM_F_OUT_SAM)
		for (i = 0; i < (int32_t)pl.mi->n_seq; ++i)
			printf("@SQ\tSN:%s\tLN:%d\n", pl.mi->seq[i].name, pl.mi->seq[i].len);

	kt_pipeline(2, worker_pipeline, &pl, 3);

	free(pl.str.s);
	mm_idx_destroy(mi);
	free(pl.rid_shift);
	for (i = 0; i < n_split_idx; ++i)
		fclose(pl.fp_parts[i]);
	free(pl.fp_parts);
	for (i = 0; i < pl.n_fp; ++i)
		mm_bseq_close(pl.fp[i]);
	free(pl.fp);
	mm_split_rm_tmp(opt->split_prefix, n_split_idx);
	return 0;
}
